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Abstract: We present a new formalism for parton shower simulation of QCD jets, which 
incorporates the following features: invariance under boosts along jet axes, improved treat- 
ment of heavy quark fragmentation, angular-ordered evolution with soft gluon coherence, 
more accurate soft gluon angular distributions, and better coverage of phase space. It is 
implemented in the new HERWIG++ event generator. 
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1. Introduction 



The parton shower approximation has become a key component of a wide range of com- 
parisons between theory and experiment in particle physics. Calculations of infrared-safe 
observables, i.e. those that are asymptotically insensitive to soft physics, can be performed 
in fixed-order perturbation theory, but the resulting final states consist of a few isolated 
partons, quite unlike the multihadron final states observed experimentally. One can at- 
tempt to identify isolated partons with hadronic jets, but then the energy flows within and 
between jets are not well represented. 

At present, the only means of connecting few-parton states with the real world is via 
parton showers, which generate high-multiplicity partonic final states in an approximation 
that retains enhanced collinear and soft contributions to all orders. Such multiparton 
states can be interfaced to a hadronization model which does not require large momentum 
transfers in order to produce a realistic hadronic final state. Hadronization and detector 
corrections to the fixed-order predictions can then be computed, and the results have 
generally been found to be in satisfactory agreement with the data. Infrared-sensitive 
quantities such as hadron spectra and multiplicities have also been described successfully 
using parton showers. This increases confidence that similar techniques can be used to 
predict new physics signals and backgrounds in future experiments. 

A crucial ingredient of modern parton showering algorithms 1 is angular ordering, which 
ensures that important aspects of soft gluon coherence are included in an azimuthally- 
averaged form. The angular shower evolution variable [2] used in the event generator 
program HERWIG [3] is good for ensuring that angular ordering is built in from the outset, 
but the phase space is complicated and not invariant under any kind of boosts. Evolution 
in virtuality looks natural but then angular ordering must be imposed afterwards, as is 
done in PYTHIA [4]. 

In the present paper we investigate a new shower evolution formalism, based on an 
angular variable related to transverse momentum [5-8] . The main aim is to retain the direct 
angular ordering of the shower while improving the Lorentz invariance of the evolution and 
simplifying the coverage of phase space, especially in the soft region. The new shower 
variables also permit a better treatment of heavy quark fragmentation, through evolution 
down to zero transverse momentum and the use of mass-dependent splitting functions, 
which eliminate the sharply-defined collinear "dead cones" seen in earlier angular-ordered 
treatments. 

In the following section we define the new shower variables and their associated kine- 
matics and dynamics, including the appropriate argument of the running coupling, the 
mass-dependent parton branching probability, and the shower evolution cutoff. The vari- 
ables are defined slightly differently for initial- and final-state parton branching, and depend 
on the colour connection of the evolving parton, so we consider in subsequent sections the 
various possible configurations of colour flow between initial and final jets. 

The formalism presented here is implemented in the new Monte Carlo event generator 
HERWIG++ [9]. Results for e + e~ annihilation and comparisons with LEP data will be 

1 For a general introduction to the parton shower approximation, see for example Chapter 5 of [1]. 
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presented shortly in a separate publication [10]. The formulae in the present paper could 
also be used to construct a matching scheme for next-to-leading order (NLO) QCD cal- 
culations and HERWIG++ parton showers, similar to that developed for HERWIG showers 
in [11,12] and implemented in the MC@NLO event generator [13], or to improve schemes 
for combining parton showers with fixed-order matrix elements [14]. 

2. New variables for parton branching 

2.1 Final-state quark branching 
2.1.1 Kinematics 

Consider parton branching in an outgoing (heavy) quark jet. Define the quark momentum 
after the ith gluon emission q^\ — > qi + ki (see fig. [l]) as 

Qi = OLip + fan + q±i (2.1) 

where p is the jet's "parent parton" momentum (p 2 = m 2 , the on-shell quark mass-squared), 
re is a lightlike "backward" 4- vector (n 2 = 0), and q±i is the transverse momentum (g^- = 
-q±i , Q±i ■ P = q±i ■ n = 0). Then 

ft = ^? + «?-*"' ■ (2.2) 
Lctip ■ n 




Figure 1: Final-state parton branching. The blob represents the hard subprocess. 

The momentum fraction and relative transverse momentum are now defined as 

Zi = , p ±i = q ±i - Ziq ±i _i . (2.3) 

OLi-l 

Then we have 

n 2 k 2 n , 2 

llx= q t + T377 + 7iTr7V (2 ' 4) 

Z{ 1 Zi Zi y 1 Zi) 

2.1.2 Running coupling 

To find the optimal argument of a s , we consider the branching of a quark of virtuality 



q 2 into an on-shell quark and an off-shell gluon of virtuality k 2 [15]. From eq. (2.4), the 
propagator denominator is 

z 1 — z z(l — z) 1 — z \ z 
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The dispersion relation for the running coupling is supposed to be 



fi 2 /i 2 J k 2 (k 2 + n 2 

where Ps(k 2 ) is the discontinuity of a s (— k 2 ). The first term on the right-hand side comes 
from cutting through the on-shell gluon, the second from cutting through the gluon self- 
energy. In our case we have k 2 + [p_|_ 2 + (1 — z) 2 m 2 ]/z in place of k 2 + /x 2 . We are interested 
in soft gluon resummation (z — > 1) [8] and so we ignore the factor of \/z here. Thus the 
suggested argument of a s is p_|_ 2 + (1 — z) 2 m 2 . In practice we impose a minimum virtuality 
on light quarks and gluons in the parton shower, and therefore the actual argument is 
slightly more complicated (see below). 

2.1.3 Evolution variable 

The evolution variable is not simply q 2 since this would ignore angular ordering. For 
massless parton branching this means the evolution variable should be p^ 2 /[z(l — z)] 2 = 
q 2 /z(l — z) [18]. For gluon emission by a massive quark we assume this generalizes to 
(q 2 — m 2 )/z{l — z). To define a resolvable emission we also need to introduce a minimum 



virtuality Q g for gluons and light quarks. Therefore from eq. (2.5) the evolution variable 



n 

is 

a 2 - P ^ 2 i ^ i Q ° (27) 
q z 2 {l-z) 2 z 2 z(l-z) 2 [ ' 

where fi = max(m, Q g ). For the argument of the running coupling we use 

z 2 (l - z) 2 q 2 = p ± 2 + (1 - z)V + zQ] • (2-8) 

Note that for massive quarks this allows us to evolve down to p^ < (1 — z)m, i.e. inside 
the dead cone [16,17]. 

Angular ordering of the branching qi — > qi + \ is defined by 

q i+ i < Ziqt . (2.9) 

The factor of Z{ enters because the angle at each branching is inversely proportional to the 
momentum fraction of the parent. Similarly for branching on the gluon, ki — > ki+\, we 
require 

k i+ l < (1 - Zi)q t . (2.10) 

2.1.4 Branching probability 

For the parton branching probability we use the mass-dependent splitting functions of 
ref. [7]. These are derived in the quasi- collinear limit, in which p_|_ 2 and m 2 are treated as 
small (compared to p ■ n) but p_|_ 2 /m 2 is not necessarily small. In this limit the q — > qg 
splitting function is 



P gq (z, P± 2 ) = C F 



1 + z 2 2z(l - z)m 2 



p_i_ 2 + (1 — z) 2 m 2 



(2.11) 
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Note that at = the factor in square brackets is just 1 — z, i.e. the soft singularity 
at z — > 1 becomes a zero in the collinear direction. The minimum virtuality Q 2 serves 
only to define a resolvable emission, and therefore we omit it when defining the branching 
probability in terms of the evolution variable (|2.7D as 



dP{q qg) = P qq dz = ^a s [z 2 (l 

2.2 Gluon splitting 



_2 dq dz 



\2~2 



n2 



1 



1 + z 2 



2m 2 

zq 2 



(2.12) 



In the case of a final-state gluon splitting into a pair of heavy quarks of mass m, the 
quasi-collinear splitting function derived in [7] is 



P q9 (z, P± 2 ) = T R 



1 - 2z(l 



p ± 2 + m 2 



(2.13) 
i at the phase 



We note that this splitting function is bounded above by its value Tr 
space boundary p_|_ = 0, and below by Tr/2. By analogy with eq. (2.7), in this case 
the evolution variable q is related to the virtuality of the gluon or the relative transverse 
momentum of the splitting by 



p_l_ 2 + m 2 



z(l-z) z 2 (l-z) 2 ' 
In terms of the variables q, z, the g — > qq branching probability then reads 

2m 2 



dP(g 



T R r 2/-, \2~2idq 

—a s [z (l-z) q ]-^- 



1 - 2z(l - z) + 



z(l 



z)q 2 



dz . 



(2.14) 



(2.15) 



In the case of gluon splitting into gluons, the branching probability takes the familiar 



form 



dP(g -> gg) 



Ca 

2ir 



a s [z 2 (l 



z) q 1 
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+ 



+ z{l-z) 



dz . 



(2.16) 



Since we introduce a minimum virtuality Q 2 for gluons, the relationship between the evo- 
lution variable and the relative transverse momentum for this splitting is as in eq. ( 2.1 j| ) 
but with the heavy quark mass m replaced by Q g . Similarly, for gluon splitting to light 
quarks we use eq. (2.14) with [i = max(m, Q g ) in place of m. 

2.3 Initial-state branching 

Consider the initial-state (spacelike) branching of a partonic constituent of an incoming 
hadron that undergoes some hard collisions subprocess such as deep inelastic lepton scat- 
tering. The momenta are defined as in eq. {2A), with the reference vector p along the beam 
direction. In this case the evolution is performed backwards from the hard sub-process to 
the incoming hadron, as shown in fig. |2[ Thus we now define in place of eq. (2.3) 



CLU-1 — z iQ±_i 



(2.17) 
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Figure 2: Initial-state parton branching. The blob represents the hard subprocess. 



Then 



Qi-i 



d - ^-k} - . (2.18) 

We assume a massless variable-flavour-number evolution scheme [19,20] for constituent 
parton branching, setting m = and putting all emitted gluons at the minimum virtuality, 
kf = Q 2 . The angular evolution variable now relates only to the angle of the emitted gluon 
and therefore we choose 

4? - P -f^f , (2.1.) 

with ordering condition simply qi + \ < qi. Correspondingly, for the argument of the running 
coupling we now use (1 — z) 2 q 2 . 

A different type of initial-state branching occurs in the decay of heavy, quasi-stable 
coloured objects like the top quark. Here the momentum of the incoming heavy object 
is fixed and evolution is performed forwards to the hard decay process. In this case we 
cannot neglect the mass of the parton and eq. ( [2.19 ) becomes 



P±i 2 + ZiQ 2 



9 +m 2 



(1 " z t ) 2 

while the branching probability ( 2.12; ) is replaced by 

dz 



Z7T q 



2 1 



l + z 2 



2zm 



21 



(2.20) 



(2.21) 



2.4 Allowed regions and termination of branching 

The allowed phase space for each branching is given by requiring a real relative transverse 
momentum, p_|_ 2 > 0. In final-state q — > qg branching, we have from eq. ( |2.8[ ) 



z 2 (l-z) 2 q 2 >(l-z) 2 l , 2 + zQ 2 g . 
This yields a rather complicated boundary in the (q, z) plane. However, since 

(1 _ z )2y + zQ 2 > (1 _ z)2/i 2 Z 2 Q 2 



(2.22) 



(2.23) 
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we see that the phase space lies inside the region 

— < z < 1 - , (2.24 
q q 

and approaches these limits for large values of q. The precise phase space can therefore 
be filled efficiently by generating values of z between these limits and rejecting those that 
violate the inequality ( 2.22| ). The resulting threshold for q is slightly larger than but of the 
order of m + Q g . 

In gluon splitting, we obtain the allowed phase space range from eq. ( 2.14| ) as 



z_ < z < z+, z± = ^ (l ± y 1 - -^j and q > 4/i (2.25) 

where fi = m for splitting into heavy quarks, or [i = max(m, Q g ) more generally. Therefore, 
analogously to eq. Q2.24 ), the phase space lies within the range 

- < z < 1 - - . (2.26) 
q q 

Schematically, the parton shower corresponds to selecting a sequence of (<ji,Zi) values 
by solving the equations 



rqi-i 




= dz-— dz— (2.27) 



where S [0, 1] are uniform pseudorandom numbers. Whenever the algorithm selects 
a value of q below the threshold, branching of that parton is terminated. The minimum 
virtuality Q g thus determines the scale at which soft or collinear parton emission becomes 
unresolvable. In the absence of such a scale one eventually reaches a region where the 
perturbative expression for the running coupling is divergent. 

One may wish to use a parametrization of a s at low scales such that a s (0) is finite. 
However, a cutoff Q g is still needed to avoid divergence of the q — > qg and g — > gg branching 
probabilities. Alternatively one could consider parametrizing a s such that a s (0) = 0, e.g. 

a S (q 2 ) = -^ra s (Ql) for q < Q c , (2.28) 
where Q c > A. Then the total branching probability below Q c is (for massless quarks) 

P c (q - qg) = C F ^^ f 1 2,(1 + z 2 ) dz = , (2.29) 

2vr J tt 

and no explicit cutoff is required, although of course Q c is essentially playing the same role. 

After branching has terminated, the outgoing partons are put on mass-shell (or given 
the virtual mass Q g if lighter) and the relative transverse momenta of the branchings in 
the shower are computed. For final-state gluon splitting we have 



|pj = V* 2 (1-*)V-M 2 , (2-30) 
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or else, if the parent is a quark, 



IpJ = ^{l-zY{z^-^)-zQl (2.31) 

The virtualities of the internal lines of the shower can now be computed backwards accord- 
ing to eq. fl2.4|). Finally, the azimuthal directions of the Pj_'s can be chosen [21] and the 



full 4- momenta reconstructed using eqs. (2.1) and (^ 



In initial-state constituent parton branching the evolution is "guided" by the parton 
distribution functions (PDFs) of the incoming parent hadron. Since PDFs are often not 
tabulated below some scale Q s > Qq, one may wish to terminate branching whenever 
q < Q s is selected. In that case the incoming parton is assigned virtuality q\ ~ — Q 2 and 
the spacelike virtualities of internal lines are then reconstructed back from q 2 to q$ using 
the transverse momenta deduced from eq. ( |2.19| ) inserted in eq. ( |2.18| ). 

For initial-state branching in the decay of a heavy, quasi-stable coloured object, the 
branching proceeds in the opposite direction but the reconstruction of momenta is similar, 
using eq. ( pT20| ) instead of ( gig ). 



2.5 Treatment of colour flows 

The more detailed treatment depend on the choice of the "backward" vector n and on 
which quantities are to be held fixed during jet evolution. Normally n should be taken 
along the colour-connected partner of the radiating parton, and the 4-momentum of the 
colour-connected system should be preserved. The upper limits on the evolution variable 
q for the colour-connected jets should be chosen so as to cover the phase space in the soft 
limit, with the best possible approximation to the correct angular distribution. In setting 
these limits we neglect the minimum virtuality Q 2 , which is a good approximation at high 
energies. We consider separately the four cases that the colour connection is between two 
final-state jets, two initial-state (beam) jets, a beam jet and a final-state jet, or a decaying 
heavy parton and a decay-product jet. 



3. Final-final colour connection 



Consider the process a — ► b + c where a is a colour singlet and b and c are colour-connected. 
Examples are e + e~ — ► qq and W — > qq'. We need to preserve the 4-momentum of a and 
therefore we work in its rest-frame, 

Pa = Q(l,O,0), p b = ±Q(l + 6-c,0,A) , p c = |Q(l-6 + c,0,-A) , (3.1) 

where pi = Q 2 , b = ml/Q 2 , c = m 2 /Q 2 and 

A = A(l, b, c ) = y/(l + b- c) 2 - 46 = 7(1 - b + c) 2 - 4c . (3.2) 

For emission of a gluon g from b we write 

qi = aiPb + Pin + q ±i (3.3) 
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where q_|_ 9 = kj_, q_|_ 6 = — kj_, q_|_ c = and we choose 

n = ±Q(A,0,-A) . (3.4) 

Notice that, if c is massive, the alignment of n along p c is exact only in a certain class 
of Lorentz frames. However, if we try to use a massive "backward" vector the kinematics 
become too complicated. 

To preserve p a = Qb + Qc + Qg we require 

Va ; = Vft = — (3.5) 

^ ^ H l + b-c + X y ' 

whereas the mass-shell conditions give 

P b = \n i u rTT ba i 

A(l + b — c + X) \ ab 

ft= A(1+i , 2 _ e+A) (| -fa.) (3-6) 

f ' s = A(l+f>-e + A) " 
where k = k± 2 /Q 2 - Our new variables are 

^ ~ Q 2 7 MM 

z = ■ , k = > b , (3.7) 



where from eq. (2.7) we have 

K = {z 2 R-b){\- z) 2 } (3.8) 
and so \JbjR < z < 1. From eqs.( |3.5|) -( |3.8D we find 

z 



l+6-c+ A 



1 + b - c + z(l - z)k + v^[i - b + c- z(l - z)R] 2 - 46 



c*c= , 2 (3-9) 
l + o — c + A z 

1 - 2 

a 9 = a 6 , 

z 

with the /3j's given by eq. (j3l|). 
3.1 Phase space variables 

It is convenient to express the phase space in terms of the Dalitz plot variables 

Xi = = (1 + b - c)oi + XPi . (3.10) 

Substituting from eqs. (|3.6| ) and (^S|), we find 

x c = 1 — b + c — z(l — z)R 

x b = (2 - x c )r +(z- r)^Jxl - 4c (3.11) 



(2 - x c )(l - r) - (z - r)^Jx 2 c - 4c 
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where 

r= H l+ i+^.'- <3 ' i2) 

The Jacobian factor is thus simply 

d(x b , x c 



z(l - z)yjx 2 c -±c (3.13) 



d(z, R) 

and the quasi-collinear branching probability ( 2.1 2| ) translates to 

a s dx b dx c \l + z 2 2b 

dP(q -> qg) = C F — - 

27r (l - b + c - x c ) \Ac - 4c 

where 



1 — z I — b + c — x 



(3.14) 



— (2 — x r )r 

z = r+ 7^== ^, (3.15) 

r being the function of x c given in eq. ( 3.12j ). 
For emission from parton c we write 

Qi = aiPc + Pin + q±i (3.16) 

where now we choose 

n = iQ(A,0,A). (3.17) 

Clearly, the region covered and the branching probability will be as for emission from 
parton b, but with Xf, and x c , b and c interchanged. 

3.2 Soft gluon region 

For emission from parton b in the soft region 1 — z = e — >0we have 

x c ~l — fe + c — en , Xfe ~ 1 + 6 — c — eh' (3.18) 

where 

k' = \+^(l-b-c-\) . (3.19) 

Since k is an angular variable, we can express it in terms of the angle 9 bg between the 
directions of the emitting parton b and the emitted gluon in the rest frame of a. In the 
soft region we find 

_ (l + b-c + X)(l + b-c- XcosObg) 

K = ; — : — (3.2U) 

2(1 + cos 9 bg ) v ' 

Thus R = b at 6 bg = and R — > oo as 9 bg — > n. 

For soft emission from parton c, the roles of x b and x c , b and c are interchanged. To 
cover the whole angular region in the soft limit, we therefore require R < R b in jet b and 
R < R c in jet c, where 

%=~^ (3.21) 

and hence 

(R b -b)(R c -c) = ^{1-b-c + X) 2 . (3.22) 
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In particular, the most symmetric choice is 



K b 



i(l + 6-c + A) 



i(l-6 + c + A) . 



(3.23) 



The largest region that can be covered by one jet corresponds to the maximal value of k 
allowed in eq. (|3.9| ) for real o^, i.e. for the maximal b jet 



Kb 



4(1 -2Vb-b + c) 



(3.24) 



3.3 Example: e + e 

Here we have b = c 

e + e~ — > 



qqg 



p, A = -^/l — 4/3 = the quark velocity in the Born process 
The phase space and the two jet regions for the symmetrical choice ( 3.23j ) are 
shown in fig. ||. The region D, corresponding to hard non-collinear gluon emission, is not 
included in either jet and must be filled using the 0(a s ) matrix element (see below). 
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Figure 3: Phase space for e + e — > qqg for m q = 5 GeV, Q 2 = 171%, with symmetric definition of 
quark and antiquark jets. 

For the maximal quark jet we get from eq. ( |3.24j ) 

k q = 4(1 - 2Jp) , (3.25) 
as shown in fig. || together with the complementary antiquark jet region given by eq. (|3,22| ). 
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Figure 4: Phase space for e + e — > qqg for m q = 5 GeV, Q 2 = m 2 z , with maximal region for the 
quark jet. 



3.3.1 Exact matrix element 

The e + e~ — > V — > ggg differential cross section, where 1/ represents a vector current such 
as a virtual photon, is given to first order in a s by [22, 23] 



1 d 2 ay a 3 Cf 



oy dx 9 dxq 2ir v 
where 

and 



2p 



2p 



(x q + 2pf + (Xg + 2p) 2 + Cy 

(l + 2p)(l-x,)(l-x ? ) (1-x,) 2 (1-x,-) 2 

Cy = -8/9(1 + 2p) 



0V = co (1 + 2p) f 



(3.26) 



(3.27) 



(3.28) 



is the Born cross section for heavy quark production by a vector current, ctq being the 
massless quark Born cross section. 

In the case of the axial current contribution e + e~ — > A — ► qqg, instead of eq. ( 3.2C| ) 
we have 



1 cPcta «s Cf 



a a dx„ dx^ 2ir v 



(x q + 2p) 2 + ( X g + 2p) 2 + ( A 

V 2 (l - Xg)(l - Xq) 



2 P 



2 P 



[1 - Xq) 2 (1 - Xq) 2 



, (3.29) 
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where 

( A = 2p[(3 + x g f - 19 + 4p] , (3.30) 
a a being the Born cross section for heavy quark production by the axial current: 

a A = (T v 3 . (3.31) 

3.3.2 Soft gluon distribution 

In the soft gluon region 1 — z = e — > the branching probability (|3.14|) becomes 

d 2 P a s 2C F 



dx g dxq 2ir ve 2 



■/.(«) 



/.(«) = 4- 4- (3.32) 



In this limit, the exact vector and axial current matrix elements, eqs. ( p.26| ) and ( |3.2 
respectively, give identical distributions: 

1 d 2 a v 1 d 2 

/(«) 



cy dxg dx^ <T^4 2ir ve 2 

l-2p p p 



v 



R 2 k' 2 

2 



/.(«) • (3-33) 



Since from eq. (|3.19| ) 



= « + - — - ) > v , (3.34) 



the parton shower approximation ( p.32| ) always overestimates the true result in the soft 
limit, and so correction by the rejection method is straightforward. For small values of p 
we have 

f(K) = l-^ + 2 -^-2p + 0{p 2 ) . (3.35) 

K K K 

Since k > p we see that the error in the approximation ( |3.32| ) is at most O(p), for any 
value of k (fig. |5|). 

3.3.3 Dead region contribution 

The integral over the dead region may be expressed as 

— / d 2 a v = ^C F F$(k q ) (3.36) 

where k q parametrizes the boundary of the quark jet. As shown in fig. ^, this is actually 
maximal, but still small, at the symmetric point given by eq. ( p. 23 ). 

Although the integral in eq. ( 3.36 ) is finite, the integrand diverges as one approaches 
the soft limit x q = Xq = 1 via the narrow "neck" of the dead region in fig. || or |||. This 
could cause problems in generating qqg configurations in the dead region in order to apply a 
matrix element correction [24]. To avoid such problems, one can map the region x q , x q > | 
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Figure 5: The function f(k) giving the gluon angular distribution in the soft limit, for m = 5 
GeV, Q 2 = m|. The exact result eq. ( 3.33 ), solid curve, and shower approximation ( 3.32 ), dashed, 
are not distinguishable on this scale. 



into a region whose width vanishes quadratically clS X q - X q — * 1, as illustrated in fig. 0. The 
mapping shown is 



4-(l-* 9 ) 



7 

4 " ' 



q > 



l-2(l-< 



7-(!-^) 



5 1 

2' 



(3.37) 



when Xq > Xq > j. Within the mapped region, the integrand then has an extra weight 
factor of 2(1 — x' ) which regularizes the soft divergence. When Xq > x q > |, x q and Xq are 
interchanged in both the mapping and the weight. 



4. Initial-initial colour connection 

Here we consider the inverse process b + c — > a where a is a colour singlet of invariant 
mass Q and b, c are beam jets. The kinematics are simple because we take beam jets to be 
massless: in the cm. frame 

p a = Q(l,O,0), p 6 = iQ(l,0,l) , p c = iQ(l,0,-l). (4.1) 
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For emission of a gluon g from b we write 

Qi = aip b + (3ip c + q ±i (4.2) 

where q_|_ 9 = kj_, q ±a = — k^, q_|_ fe = q_|_ c = 0. Notice that in this case the recoil transverse 
momentum is taken by the colour singlet a so we cannot preserve its 4-momentum. We 
choose to preserve its mass and rapidity, so that 



a a = p a = VI + « , (4.3) 
where as before k = k± 2 /Q 2 . Now we have 

Pb = a c = , a g P g = K , 

a a = a b - a g , (3 a = /3 C - P g , (4.4) 
and our new variables in this case are 

1 a 9 ~ ? 2 K I A X\ 

a b Q 2 (1 - z) 2 
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Figure 7: The soft region, with jet boundaries (solid) and mapped region (dashed), for m = 5 
GeV, Q 2 = m%. 

Thus we find 

a a = Pa = y/l + (l- z) 2 k , 



a b = -y/l + {l-z) 2 k, 

4.1 Phase space variables 

It is convenient to express the kinematics in terms of the "reduced" Mandelstam invariants: 

s = (q b + q c ?/Q 2 , t = (q b -q g ) 2 /Q 2 , u = (q c - q g f /Q 2 . (4.7) 

The phase space limits are 

l<s<SyQ 2 , l-s<i<0, u=l-s-i (4.8) 

where S is the beam-beam cm. energy squared. In terms of the shower variables for beam 
jet b, we have 

s = a b p c = -[1 + (1 - z)R] , i = -a b p g = -(1 - z)k , u = -(l-z)s. (4.9) 
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Figure 8: Beam jets (B,C) and dead region (D) in initial-state branching. 

Thus curves of constant k in the (s, t) plane are given by 

_ k(l-s) 



(4.10) 



K + S 

and the Jacobian factor for conversion of the shower variables to the Mandelstam invariants 
is 



Ir4 = — *■ < 411 > 

o(z, K) z 
For the other beam jet c we have t <-> u and thus 

f= g(1 " g) . (4.12) 

K + S V ' 

We see that in order for the jet regions to touch without overlapping in the soft limit s — > 1, 
t — > 0, we need k < Kb in jet b and k < k c in jet c, where k c = 1/kf,. The most symmetrical 
choice = 1, as shown in fig. || but we can take kb or k c as large as we like. 

4.2 Example: Drell-Yan process 

Consider radiation from the quark in the Drell-Yan process, qq — > In the laboratory 
frame we have 

q q = (Px h 0, Px x ) , qg = {Px 2 , 0, -Px 2 ) (4.13) 

where P = is the beam momentum. If we generated the initial hard process qq — > 
with momentum fractions and we want to preserve the mass and rapidity of the Z° 

we require 

x\ = x q a b , x 2 = Xqf3 c (4.14) 

where and (5 C are given by eqs 
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The branching probability in the parton shower approximation is 



* P ft*^, (4.15) 



dz dk 2tt k 1 

which gives a differential cross section (s = sQ 2 , etc.) 

1 d 2 a D(xi )D(xo) «s „ s + u r9 N9n . 

r~r= in) ( tt-Cf^— \s 2 + (s + u) 2 ] (4.16) 

a dsdt D(x q )D(xg) 2vr * s 3 tu L v 7J v ' 

where ctq is the Born cross section. The functions D(x\) etc. are parton distribution 
functions in the incoming hadrons; these factors take account of the change of kinematics 
x q ,xq — > xi,X2 discussed above. 

The exact differential cross section for qq — > gZ° to order a s is 

^^ = ^(x 9 )^-) 2^ g ^[ (g + t) + • 



Since Q 2 = s + t + u and t < 0, we see that the parton shower approximation ( 4.16j ) 



overestimates the exact expression, becoming exact in the collinear or soft limit t —> 0. 
Therefore the gluon distribution in the jet regions can be corrected efficiently by the re- 
jection method, and the dead region can be filled using the matrix element, as was done 
in [25]. The benefit of the new variables is that the angular distribution of soft gluon 
emission requires no correction, provided the jet regions touch without overlapping in the 
soft region. As shown above, this will be the case if the upper limits on k satisfy kq = l/k q . 



5. Initial-final colour connection 



Consider the process a + b — > c where a is a colour singlet and the beam parton b and 
outgoing parton c are colour-connected. An example is deep inelastic scattering, where a 
is a (charged or neutral) virtual gauge boson. We need to preserve the 4-momentum of a 
and therefore we work in the Breit frame: 

Pa = Q(0,0,-1) , p 6 = iQ(l + c,0,l + c) , p c = ±Q(l + c,0,-l + c) , (5.1) 

where p 2 = —Q 2 , p\ = 0, and m 2 = cQ 2 . Notice that the beam parton b is always taken 
to be massless, but the outgoing parton c can be massive (e.g. in W + d — ► c). 

5.1 Initial-state branching 

For emission of a gluon g from the incoming parton b we write 

qi = otip h + fan + qxi (5.2) 
where q_|_ 9 = k^, q ±b = 0, q_|_ c = — k^ and we choose 

n = iQ(l + c,0,-l -c) . (5.3) 
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To preserve p a = q c + q g — q b we now require 

a b - a c - a g = 0c + (5 g - (5 b 
whereas the mass-shell condition is 

a 4 /3 4 Q 2 (l+c) 2 =q ±i 2 + g 4 2 

which gives 

c (I 1 

1 + C= — + K — + — 



1 + c 



(5.4) 



(5.5) 



a r a, 



The new variables for emission from the beam jet are as in eq. 
5.6T), we find 



(5.6) 

Substituting in 



1 



"6 



2z{l + c) 



a c = za b 



1 



1 + c + (1 - z)k + + c + (1 - z)k] 2 - - z)k 

, a g = (1 - z)a 6 , ft = , 
(1 - z)k 



1 + c 

1 C+(1-Z) 2 K 



l + c' z{\ + c)a b - 1 ' ^ 9 
5.2 Final-state branching 

Next consider emission from the outgoing parton c. In this case we write 

Qi = aip c + f3ip b + q ±i 
To preserve p a = q c + q g — q b we require 

a c + a g ~ oc b = ft - ft - P g = 1 
whereas the mass-shell condition is now 

otiPi {Q 2 + m 2 ) = q ±i 2 + q 2 - a 2 m 2 c . 



(1 + cfa b 



Or 



a c + a c 



<r_]_ 

Q 2 ~ 'z 2 



c + 



{l-zf 



Thus in this case we find 



a b = , a c = z , a g = 1 — z , 
1 



ft 



l + c-M(l- z)k , 



1 + c 1 
1 — z 

ft = -r— [2c + z(l - z)k] 



l + c l 

1 ~ z r^,. 
1 + c 



z l k - 2c] . 



(5.7) 



(5.8) 



(5.9) 



(5.10) 



The new variables for emission from an outgoing parton are as in eqs. (3.7,3.8) with b 
replaced by c: 



(5.11) 



(5.12) 
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5.3 Phase space variables 

In this process the invariant phase space variables are usually taken to be 

_ Q 2 _ 1c ■ Qb 

X P — r) ) Z P ~ 

2p a -q b ^ Pa-Qb 
In terms of the new variables for emission from the beam parton, we have 



(5.13) 



1 



(1 + c)a b 



2z 1 + c + (1 - z)k + y/[l + c + (1 - z)kf - 4z(l - z)k) (5.14) 



(1 + c)p c = - (l - c - (1 - z)k + V[l + c + (1 - z)k] 2 - 4z(l - zjkj , (5.15) 



with the Jacobian 



d(xp,z p ) = lfj_ + 1 + c _ 2 

d(z,k) 



(5.16) 



In the soft limit z = 1 — e we therefore find for the beam jet 

1 ecK 



1 + c 



1 -e 



(1 + c); 



, z p ~ 1 



€K 
T+C 



and 



d{xp, Zp) 



d(z,K) (1 + c) 2 ' 

In terms of the variables for emission from the outgoing parton, 

1 1 



(5.17) 



(5.18) 



(l + c)/? 6 1 + c + z(1-z)k 



so the Jacobian is simply 



and in the soft limit 



d(x p , Zp) _ /, % 2 



l + c 



T+c 



, z p ~ 1 - e 



(5.19) 



(5.20) 



(5.21) 



with the Jacobian again given by eq. (gig) . For full coverage of phase space in the soft 
limit we require k < k b in jet 6 and k < k c in jet c, where 



k 6 (k c - c) = (1 + c) 2 . 



(5.22) 



Thus the most symmetrical choice is k b = 1 + c, k c = 1 + 2c, as shown in fig. ^. On the 
other hand, any larger or smaller combination satisfying eq. ( |5,22| ) is allowed, as illustrated 
in fig. 1C for k b = 10. 
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Figure 9: Beam jet (B), outgoing jet (C) and dead region (D) in initial-final state branching: 
c = 0.25, k b = 1.25, k c = 1.5. 




0.8 



Figure 10: Beam jet (B), outgoing jet (C) and dead region (D) in initial-final state branching: 
c = 0.25, k b = 10, k c = 0.40625. 



5.4 Example: deep inelastic scattering 

Consider deep inelastic scattering on a hadron of momentum by exchange of a virtual 
photon of momentum q^ 1 . If the contribution to the Born cross section from scattering on 
a quark of momentum fraction xb = Q 2 /2P • q is represented by o"o (a function of xb and 
Q 2 ), then the correction due to single gluon emission is given by 



1 d 2 a C_pa s D(x B /x p ) 1 + (x p + z p - 1) 



a^dxpdzp 2ir D(xb) x p (1 — x p )(l — z p ) 



(5.23) 
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In the soft limit i p ,z p ->lwe have, from eqs. ( |5.17 , 5T2~l| ) with c = 0, 

(1 - x p )(l - z p ) ~ e 2 k 

and so 



1 d l o 



Cpa s 1 



whereas the parton shower approximation gives 

1 d 2 a Cpds 1 



ctq dz dk 



7T en 



(5.24) 
(5.25) 

(5.26) 



Since the Jacobian factor ( 5.16j ) or (5.20) in this limit is simply e, the shower approximation 
is exact in the soft limit. 

5.5 Example: qq — > ti 

We denote the momenta in this process by p a + pi, — > p c + pd and the 2^2 invariants by 



2pa-Pb , t = -2p a ■ p c 



u 



-2Pa ■ Pd , 



(5.27) 



so that s + t + u = 0. Colour flows from q to t and anticolour from q to t. Therefore the 
momentum transfer q = p a — p c = p^ — p^ is carried by a colour singlet and we preserve 
this 4-momentum during showering. 

For emission from the incoming light quark or the outgoing top quark, we work in the 
Breit frame for this system, where 



g = Q(0,0,l), Pa = iQ(l + c,0,l + c) 



Pc 



§Q(l + c,0,-l + c) 



(5.28) 



with Q 2 



-t — mf and c 



m: 



''/Q 2 . Then the treatment of sects. 5.1 and |5.2| can be 



applied directly, with the substitution b — > a since the emitting system is now (a, c) rather 
than (6, c). However, the phase space variables are no longer those of sect. |5^ since they 
involve the momenta of the q and i, which in the frame ( |5.2S] ) take the general form 



p b = [IQV(1 + c) 2 + 4K, Q ± , -|Q(1 + c)] , 
Pd = [iQV(l + c) 2 + 4i^, Q ± , \Q{1 - c)] , 



where K = Qj_ /Q is related to the 2^2 invariants: 



ln2 



1 + 4/1 + 



(1 + cf 



-Q 2 (l + c) 



u 



and so 



Q 



i + 



(1 + c)i 



For emission from the incoming light quark we define as in sect. 5.1 

Qi = OLip a + fan + q±i 



(5.29) 



(5.30) 



(5.31) 



(5.32) 
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for i = a,c,g, where q_|_ a = 0, q^ c = — kj_, q_|_ g = k^, and n is as in eq. (5J5). Then 
the aj's and /3j's are given by eqs. ( |5.7| ) with the substitution b — > a. The light antiquark 
and the antitop are not affected and therefore qb = Pb, Qd = Pd- This allows the complete 
kinematics of the 2^3 process to be reconstructed. The 2 — > 3 invariants can be defined 
as in ref. [12]: 

s = 2q a ■ q b , ti = -2q a ■ q c , t 2 = -2q b ■ q d , U\ = -2q a ■ q d , u 2 = -2q b ■ q c . (5.33) 
It is convenient to express n = p c — cq so that (for i = a,c,g) 

qi = {oti - cf3i)p a + (1 + c)(3ip c + q ±i . (5.34) 

Then we find 

s = a a s , h = a a p c (l + c)i , t 2 = i , ui = a a u , u 2 = f3 c (u - ct) - a c s - 2k ± ■ Q ± . 

(5.35) 

For emission from the outgoing top we use the results of sect. 5.2, again with the 
substitution b — ► a. Thus we now have for i = a,c,g 

qi = otip c + f3ip a + q±i (5.36) 

where the ctj's and fy's are given by eqs. (5.12) with b — >• a, and we find that 

s = (5 a s , ti = a c (i a i , t 2 = i , u\ = [3 a u , u 2 = a c u - [3 c s - 2k ± • Q ± . (5.37) 

Similar formulae to eqs. fl5.35|) and ( |5.37|) , with the replacements a — > b and c — > d, 
will hold for the case of gluon emission from the colour-connected (qi) system. Using 
these relations, one can study the distribution of gluon radiation in the parton shower 
approximation and compare it with the exact qq — > tig matrix element. Agreement will 
be good in the soft and/or collinear regions but there will be regions of hard, wide-angle 
gluon emission in which matrix element corrections should be applied. Alternatively, the 
above equations can be used to formulate a modified subtraction scheme for combining 
fixed-order and parton-shower results, as was done in ref. [12] for a different parton-shower 
algorithm. 



6. Decay colour connection 

Consider the process b — ► ca where a is a colour singlet and the decaying parton b and 
outgoing parton c are colour-connected. Examples are bottom quark decay, b — > cW*, and 
top decay, t — * bW. Here we have to preserve the 4- momentum of the decaying parton b 
and therefore we work in its rest frame, 

p b = m b (l,0,0) , p c = \m b (l - a + c, 0,A) , p a = \m b (l + a - c, 0, -A) , (6.1) 

where a = m^/m?, c = m^/m^ and now 

A = A(l, a, c) = y^l + a - c) 2 - 4a = ^(1 -a + c) 2 -Ac. (6.2) 
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6.1 Initial-state branching 

For emission of a gluon g from the decaying parton b we write 



q% = aip b + fan + q±i (6.3) 

where q_|_ 9 = k^, q ±c = — k^, q ±b = and we choose 

n = §m 6 (l,0,l) , (6.4) 

i.e. aligned along p c in the rest frame of b. The mass-shell conditions give 

a a a C + K a K (ar.\ 

Pa = a a , p c = a c , p g = ctg , (6.5) 

Ol a Oi C OLg 

with k = k^ 2 /m?. From momentum conservation 

a c + K K 

a a + a c + a g = 1 1 = 1 . (6.6) 

OL a Ot C OLg 

Recall that in initial-state branching of a heavy object our new evolution variable is given 
by eq. ( |2.20| ), so we have 

a g = l-z, k = (k - 1)(1 - zf (6.7) 
where R = q 2 /ml > 1. Introducing for brevity the notation 



w = 1 — (1 — z)(k — 1) , u = I + a — c— (1 — z)k , v = \/ u 2 — Aawz , (6-8) 
from eq. ( |6.6| ) we find 

u + v i (p. n\ 

o-a = — , ol c = 1 - a a - OLg = z - a a . (6.9) 

2w 

6.2 Final-state branching 

For radiation from the outgoing parton c we write 

qi = aip c + Pin + q ±i (6.10) 

where p c is given by eq. (|6.1| ) . Since the colour-connected parton b is at rest in our working 
frame of reference, the choice of the light-like vector n in this case is somewhat arbitrary. 
By analogy with the cases treated earlier, we choose it to be opposite to that used for the 
radiation from b, i.e. along the direction of the colour singlet a: 

n = |m 6 (A, 0, -A) . (6.11) 

The kinematics are then identical with those for final-final connection (sect. ||), with the 
replacement b — > c, c — > a. 
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6.3 Phase space variables 

As in sect. ||, it is convenient to use the Dalitz plot variables, which in this case are 

2g* • Pb 



1 ~ i 



(6.12) 



For emission from the decaying parton b we have Xi = 2a.i + and hence, from eq. 
u + v u — V 



3, n 



+ 



2w 2z 
with the Jacobian factor 

<9(X a , Xg) 



, x c = w + z — x a , x g = 2 — w — z = (1 — z)R , (6.13) 



(1-z) 



d(z, R) 

In the soft limit z — > 1 — e we find 



u + v u — v a(w — z) 



21 



2w 2 2z 2 



+ 



vwz 



where 



x a ~l + a — c — eR b , x g ~ €Kb 



4 = A+^(l-a + c-A) 



(6.14) 



(6.15) 



(6.16) 



For emission from the outgoing parton c we have, from eq. (3.11) with the replacement 
b — ► c, c — > a: 



where 



x a = 1 + a — c — z(l — z)k 

x c = (2 - x a )r + (z - r)^x 2 a - 4a 

x 9 = (2 - z a )(l - r) - (z - r)^fx 2 a - 4a 



r = - 1 + 



1 + a — x c 



In the soft limit we have from eq. ( 3.1S| ) 



where 



x a ~l + a — c — ek c , x„ ~ ek c 



R' = A + — (1 - a + c- A) 
2c 



For full coverage of the soft region we require 



which gives in this case 



Kb K c 

R'b Rc 



(R b - !)(««. - c) = -(1 - a + c + A) 2 



(6.17) 

(6.18) 

(6.19) 
(6.20) 

(6.21) 
(6.22) 



Note that, while there is no upper limit on Rb, the largest value that can be chosen for R c 
is given by the equivalent of eq. ( 3.24[) , 



R c < 4(1 + a - 2yfc - c) . 



(6.23) 
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6.4 Example: top decay 

In the decay t — > Wbg we have a = m^/mf = 0.213 and c = m 2 /m 2 = 0.026, so for 
simplicity we neglect c. Then for radiation from the top we have from eq. ( |6.13| ) 

c g = (l- z)k , 



x w 



U + V u — V 



2w 2z ! 

where u, v, w are given by eqs. fl6.8|) with c —> 0. The phase space is the region 



< x g < 1 — a 



1 



x g + 1 

1 — x, 



< xw < 1 + a 



(6.24) 



(6.25) 



Notice that for real x\y we require u 2 > 4awz, i.e 

1 < k < 1 + a 1 

and 



1 



1 -a 



z(l - a) 
a{\ - z) 



< z < 1 . 



(6.26) 



(6.27) 



k + 2 v / a(K - 1) 

Thus there is no upper limit on k, but the range of z becomes more limited as k increases. 
For radiation from the b we have from eq. (|6.17| ) 



x\y = 1 + a — z(l — z)k 
x 9 = l( 2 ~ x w) ~(z-h] 



'W 



4a . 



(6.28) 



To cover the soft region we require k < kt for emission from the top quark and k < kb for 
that from the bottom, where eq. (6.22) gives 

(I -a? 



K b 



(6.29) 



kt - 1 

The most symmetrical choice would therefore appear to be kb = k t — 1 = 1 — a = 0.787, as 
illustrated in fig. 11. 

As mentioned above, there is no upper limit on k t . Thus the region covered by gluon 
emission from the top quark can be as large as we like. However, eq. ( 6.23j ) tells us that 
the upper limit for radiation from the b is 



k b < 4(1 - v 7 ^) 2 = 1.16 , 



3.30) 



and correspondingly kt > 1 + t(1 + y/o) 2 = 1.53. Figure 12 shows this maximal region 
that can be covered by emission from the 6, together with the complementary regions of 
emission from the t. 



We note from figs. |ll| and |12| that, for any value of kt, the region for emission from the 
top quark consists of two distinct parts that touch at the point x g = 1 — y/a, xw = 2y/a, 
where the W boson is at rest: a subregion Tl which includes the soft limit x g — > and a 
hard gluon region T2. 
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Figure 11: Phase space for decay t — ► Wbg, with symmetric choice of emission regions for the b 
(B) and the t (T1,T2), and the dead region (D). 



The exact t — > Wbg differential decay rate to first order in a s is given in [26]: 
1 d 2 T __ a s C F f (1 + a - x w ){\ - x g ) + x 2 g 



Tq dxwdxg 7r (1 + a — xw)x"l | 9 



1-a 



+x, 



(xw + Xg-l-a) 2 2a(l + a-x w ) 



2(1 -a) 2 



+ 



(1 -a) 2 (l + 2a) 



(6.31) 



where Tq is the lowest-order decay rate. In the soft region xw — * 1 + «> x g — > this 
becomes 



1 d 2 T a s C F 
Tq dxwdxg IT Xg 



1 



1 



(6.32) 

1 + a — xw (1 — a)x g 
For soft gluon emission (1—z = e — ► 0) from the top quark we have from eqs. fl6.15 , 6~T(f ) 



xw ~ 1 + a — e(l — a) , x g ~ ek 



(6.33) 



and so the exact form of the soft gluon distribution is 



1.3 



1.2 



1.1 



1.0 



B 



T2 
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0.0 
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Figure 12: Phase space for decay t — > VF&g, with maximal region (B) for emission from the b, 
together with complementary regions of emission from the t (T1,T2) and the dead region (D). 



where 



/*(«) 



(1 - a)k 2 ' 

In the same region the parton shower approximation ( |2.12 ) gives 



d 2 P 

dx„dxw 



1 



d 2 p 



C F 



! 1 

'1 — a)e dz dk tt (1 — a)e 2 K V « 



2"/( K ) 



(6.35) 



(6.36) 



Thus we see that, for emission from the top quark, the shower approximation is exact in the 
soft limit. At higher gluon energies, inside the region Tl the parton shower overestimates 
the exact matrix element and can therefore be corrected easily by the rejection method. 
In the hard gluon region T2, which contributes only a small finite correction to the cross 
section, the parton shower overestimates the matrix element at lower values of x g but 
underestimates it at the highest values. Therefore a combination of rejection and matrix 
element correction is needed in this region. 

For emission from the bottom quark in the soft limit, we use the results of sect. 3.2 
with the substitution 6^0, c — ► a to obtain 



xw 



1 + a 



€K 



1- a + 



(6.37) 
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Therefore the exact soft gluon distribution in the b jet should be 



1 d 2 T a s C Ff .. 

o-fb{K) 



T dxgdxw 



where 



/ 6 («) 



1 



(1 — a)k 



1 + 



(1-af 



On the other hand the parton shower approximation in this case gives simply 



d 2 P 



1 



d 2 P a, 



C F 



dxgdxw (l — a)edzdk tt (l — a)e 2 k 

Thus the soft gluon distribution in the b jet region is overestimated by a factor of 

n 2 

1 + 



(l-o)= 



(6.38) 



(6.39) 



(6.40) 



(6.41) 



which can be corrected by the rejection method. This factor varies from 1 to 5.2 for the 
symmetric choice of the b jet region kf, = 1 — a depicted in fig. 11. For the maximal b 
jet shown in fig. it rises to 8.3. Since the shower approximation is exact in the soft 
limit for emission from the top, one can reduce the amount of soft correction required 
by decreasing the b jet region and increasing that for top emission, in accordance with 
eq. Q6.29p . However, for large values of kt the dead region moves near to the collinear 
singularity at xyy = 1 + a and a large hard matrix element correction becomes necessary. 



7. Conclusions 

We have presented a new formulation of the parton-shower approximation to QCD matrix 
elements, which offers a number of advantages over previous ones. Direct angular ordering 
of the shower ensures a good emulation of important QCD coherence effects, while the 
connection between the shower variables and the Sudakov-like representation of momenta 



(2.1) simplifies the kinematics and their relation to phase space invariants. The use of 
mass-dependent splitting functions with the new variables allows an accurate description 
of soft gluon emission from heavy quarks over a wide angular region, including the collinear 
direction. The separation of showering into contributions from pairs of colour-connected 
hard partons permits a general treatment of coherence effects, which should be reliable at 
least to leading order in the number of colours. Since the formulation is slightly different 
for initial- and final-state showering, we have provided formulae for all colour-connected 
combinations of incoming and outgoing partons. 

As mentioned in the Introduction, this new shower formulation is a key element of 
the event generator HERWIG++ [9], and detailed results of its implementation for e + e~ 
annihilation will be presented shortly [10]. 
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